OVERVIEW OF MAIN STEPS IN RNA-Seq DATA ANALYSIS
(da inserire tabella)
RAW READS MANIPULATION
Accendere il computer e aprire il terminale (ctrl + alt + T)
Spostarsi nella cartella Scaricati, create una nuova cartella chiamandola tc
cd Scaricati/
mkdir tcEntrare nella cartella appena creata e valutarne il contenuto della cartella digitare il seguente comando
cd tc/
ls -lrthLa cartella dovrebbe essere vuota. Come prima cosa dobbiamo recuperare i file grezzi (FASTQ) del nostro esperimento. Per recuperare i file ci connettiamo al nostro server di laboratorio e con il protocollo di trasferimento ssh copiamo in locale (la cartella dove siamo in questo momento) i 2 file corrispondenti alle reads forward e alle reads reverse. Lavoreremo con un solo campione, ma come sapete l’esperimento consiste di 3 repliche biologiche per i non trattati e 3 repliche biologiche per quelli trattati. I comandi che qui eseguiremo possono essere ripetuti su più file contemporaneamente tramite appositi comandi come i cicli for
# prima scarichiamo il primo file (attenzione al punto in fondo al comando, indica current directory)
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_1_subset.fastq.gz .
# inserire la password quando ce la chiede ()
# controlliamo che sia stato scaricato
ls -lrth
# ora scarichiamo il secondo file
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_2_subset.fastq.gz .
# inserire nuovamente la password
# controlliamo
ls -lrthCome potete osservare la dimensione di questi file è molto bassa rispetto alla dimensione tipica di un file che risulta da un esperimento di sequenziamento come lo abbiamo descritto. Questo perché i file originali sono stati copiati e ridotti in dimensioni con lo scopo di riuscire a manipolarli e lavorarci nei tempi della lezione, principalmente per il fatto che le capacità dei computer che stiamo usando potrebbero essere limitate e quindi non in grado di eleaborare questi file in tempi brevi
Iniziamo con esplorare questi file
# "apriamo" il file
less -S Sample_Untreated_1_1_subset.fastq.gz
# contiamo il numero di reads
zcat Sample_Untreated_1_1_subset.fastq.gz | echo $((`wc -l`/4))
zcat Sample_Untreated_1_2_subset.fastq.gz | echo $((`wc -l`/4))
# il numero di reads F e R coincide?QUALITY CONTROL
Per effettuare il controllo di qualità sulle nostre reads utilizzeremo il programma FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Questo software (scritto in Java) prende come input un numero a piacimento di reads e calcola varie statistiche che possono dare un’idea della qualità dei nostri dati
# creiamo una nuova cartella dove scaricheremo il programma e salveremo i risultati
mkdir fastqc
cd fastqc/
# scarichiamo il programma con il comando wget
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
# scompattiamo la cartella
unzip fastqc_v0.11.7.zip
# rimuoviamo cartella compressa
rm fastqc_v0.11.7.zip
# vediamo cosa c'è dentro la cartella senza spostarci dentro (percorso relativo alla nostra posizione)
ls -lrth FastQC/
# il file che ci interessa, cioè il programma che andremo a richiamare è "fastqc"
# dobbiamo prima renderlo eseguibile
chmod 755 FastQC/fastqc
# controlliamo che funzioni
./FastQC/fastcq -hA questo punto possiamo lanciare il programma dando come input i due file con le reads F e R. Il comando è molto semplice, basta specificare i nomi dei file che vogliamo controllare, e la cartella dove vogliamo che vengano salvati i report coi risultati (la cartella deve già esistere quindi prima la creiamo)
# creiamo cartella per salvare i risultati
mkdir results
# ci spostiamo una cartella su, dove stanno i nostri file
cd ..
# eseguiamo il programma (dobbiamo specificare il percorso dell'eseguibile rispetto a dove siamo) e specifichiamo inoltre il percorso della cartella di output
./fastqc/FastQC/fastqc Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz -o fastqc/results/Per ogni file di input vengono creati 2 file di output, un file html e una cartella compressa. Andiamo ad esplorare il report in formato html aprendolo con Firefox
FastQC report
A questo punto commentiamo insieme i risultati del controllo qualità confrontandoci con qualità di riferimento sul sito di FastQC e con la spiegazione dei moduli di analisi (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)
Dopo aver deciso insieme quali parametri utilizzare per correggere la qualità delle nostre reads ed aver identificato gli adattatori da rimuovere, andiamo ad effettuare la rimozione degli adattatori e il quality trimming
Adapters removal and quality trimming
Per effettuare queste due operazioni useremo un unico programma che si chiama Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Il suo utilizzo è molto simile a quello di FastQC (anche Trimmomatic è scritto in Java) e necessita come input le sequenze grezze. Per ogni file di input (nel nostro caso 2 file) verranno creati 2 file di output, uno indicato come paired, l’altro come unpaired; è necessario indicare il nome di ciascuno dei file di output. Le opzioni di Trimmomatic sono:
PEindica che lavoriamo su file paired-end-trimlogindica che vogliamo salvare un file di log-threadsindica il numero di processori che vogliamo utilizzare per lanciare il programmaILLUMINACLIP:ALL_ADAPTERS.fa:2:30:10è la parte che indica la rimozione degli adattatori. Il file ALL_ADAPTERS.fa è il file in formato fasta che contiene la sequenza degli adattatori maggiormente utilizzatiLEADING:16indica la rimozione di N basi in testa alla sequenza, sotto la qualità indicataTRAILING:16indica la rimozione di N basi in coda alla sequenza, sotto la qualità indicataSLIDINGWINDOW:4:25effettua un trimming basato su una “finestra”" di basi, rimuovendo le basi nel momento in cui la qualità media per quella finestra scende sotto la qualità indicataMINLEN:36indica la rimozione delle reads che sono più corte della lunghezza indicata
# ci spostiamo nella cartella tc, creiamo nuova cartella e scarichiamo il programma
# usiamo percorso assoluto da qualsiasi punto del computer per spostarci dentro una qualsiasi cartella
cd /home/usr/Scaricati/tc # usr è specifico per il vostro utente (il mio è dlfptr70)
mkdir Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip -O Trimmomatic/Trimmomatic-0.36.zip # con l'opzione -O specifichiamo dove vogliamo scaricare il file e con che nome
# decomprimiamo senza cambiare cartella e rimuoviamo l'archivio
unzip Trimmomatic/Trimmomatic-0.36.zip -d Trimmomatic/ && rm Trimmomatic/Trimmomatic-0.36.zip
# controlliamo il funzionamento del programma
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h # -h è l'opzione per helpA questo punto possiamo lanciare Trimmomatic sui nostri file; anche se dal controllo qualità non abbiamo identificato sequenze di adattatori come contaminanti, un classico approccio potrebbe essere quello di dare come input di file degli adattatori il file fornito dal programma che quindi rimuoverà tutte le sequenze corrispondenti a quelle che trova
# creiamo file unendo in uno tutti gli adattatori
cat Trimmomatic/Trimmomatic-0.36/adapters/*.fa > Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa
# cat è un comando che unisce file "row-wise"; l'asterisco (*.fa) indica di selezionare tutti i file che hanno estensione .fa
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 10 -trimlog \
Trimmomatic/log_trimmomatic \
Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz \
Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_1_subset_trim_unpaired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_2_subset_trim_unpaired.fastq.gz \
ILLUMINACLIP:Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa:2:30:10 \
LEADING:16 TRAILING:16 SLIDINGWINDOW:4:25 MINLEN:36 &> Trimmomatic/log_trimmomatic_screen
# con il simbolo &> andiamo a salvare su un file anche l'output che normalmente sarebbe stampato sullo schermoDovrebbe terminare dopo qualche secondo senza aver stampato a schermo alcuna informazione. Andiamo a vedere se ha creato i file ed esaminiamo i file di log (log_trimmomatic e log_trimmomatic_screen); il primo contiene le seguenti info per ciascuna read:
- the read name
- the surviving sequence length
- the location of the first surviving base, aka. the amount trimmed from the start
- the location of the last surviving base in the original read
- the amount trimmed from the end
Il secondo contiene le statistiche generali
# guardiamo prima le statistiche generali
less -S Trimmomatic/log_trimmomatic_screen
# poi l'altro file
less -S Trimmomatic/log_trimmomaticA questo punto ripetiamo il controllo di qualità, questa volta sulle reads trimmate, solo sui file paired
# creiamo cartella dove salveremo nuovo controllo qualità
mkdir fastqc/results_after_trimming
# lanciamo FastQC sulle sequenze trimmate
./fastqc/FastQC/fastqc Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz -o fastqc/results_after_trimming/Andiamo ad esplorare i report e vedere se le statistiche sulla qualità sono cambiate
MAPPING (ALLINEAMENTO)
Ora i nostri file sono pronti per essere allineati sul genoma di riferimento (Arabidopsis). Per quasi tutti i programmi che effettuano allineamento il primo passaggio consiste della creazione dell’indice del genoma. Perché viene creato l’indice del genoma: citazione presa da qualcuno ma non ricordo chi (modificata) “Se entro in una biblioteca con la frase “La genetica è il ramo della biologia che si occupa del materiale ereditario” e voglio trovare il libro in cui è scritta, farò molta fatica a trovarlo, a meno che all’ingresso della biblioteca non si trovi un catalogo che mi indirizza almeno verso la sezione “Scienza””. L’indice quindi consiste in una serie di file che riassumono le principali caratteristiche del genoma di un organismo. Per effettuare l’allineamento useremo un programma chiamato STAR (https://github.com/alexdobin/STAR). Andiamo a scaricarlo insieme ai file aggiornati del genoma e del trascrittoma di Arabidopsis
# siamo nella cartella tc
mkdir STAR
# scarichiamo e salviamo nella cartella appena creata
wget https://github.com/alexdobin/STAR/archive/2.6.0c.tar.gz -O STAR/2.6.0c.tar.gz
# decomprimiamo e rimuoviamo archivio
tar -xf STAR/2.6.0c.tar.gz -C STAR/ && rm 2.6.0c.tar.gz
# controlliamo il funzionamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR -hPer creare l’indice del genoma, abbiamo bisogno del genoma. Andiamo a scaricarlo dal sito di Ensembl Plants
# genoma
mkdir genome
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -O genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# trascrittoma
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.39.gtf.gz -O Arabidopsis_thaliana.TAIR10.39.gtf.gz
# i file devono essere decompressi
gunzip genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip genome/Arabidopsis_thaliana.TAIR10.39.gtf.gzA questo punto possiamo creare l’indice con il comando di STAR --runMode genomeGenerate. Vediamo nel dettaglio le altre opzioni:
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STARchiamiamo il programma--runThreadNnumero di threads con cui vogliamo lanciare il programma--runModeindica che vogliamo creare l’indice--genomeDirindica la cartella dove vogliamo che venga salvato l’output--genomeFastaFilescartella dove è presente il file del genoma--sjdbGTFfileindichiamo che vogliamo usare anche il trascrittoma per aiutarlo nella creazione dell’indice--sjdbOverhangspecifica la lunghezza delle reads meno 1. Serve per creare il database delle giunzioni di splicing
Lanciamo il comando con le opzioni e i percorsi adattati al nostro caso
# creiamo cartella per l'indice all'interno della cartella del genoma
mkdir -p genome/index
# lanciamo la creazione dell'indice
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--runMode genomeGenerate \
--genomeDir genome/index \
--genomeFastaFiles genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--sjdbGTFfile genome/Arabidopsis_thaliana.TAIR10.39.gtf \
--sjdbOverhang 100Con 10 core dovrebbe impiegare circa 3 minuti. Andiamo a vedere se ha creato i file nella cartella index
# spostiamo il file Log.out dove dovrebbe stare
mv Log.out genome/index/
# guardiamo i file appena creati
ls -lrht/genome/index
# raccolgono diverse statistiche sul genoma di Arabidopsis, comprese info sul db di splicejunctionsPossiamo ora andare ad effettuare l’allineamento. Ovviamente per l’allineamento useremo le reads trimmate e pulite (le paired). Vediamo nel dettaglio i comandi:
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STARchiamiamo il programma--runThreadNnumero di threads con cui vogliamo lanciare il programma--genomeDirdove sta la cartella con l’indice--readFilesInqui scriviamo il percorso e il nome dei campioni da mappare--readFilesCommandqual è il formato di comrpessione dei file da mappare--outFileNamePrefixpercorso e prefisso che vogliamo assegnare ai file di output--alignIntronMaxlunghezza massima (stimata) degli introni del nostro organismo--outSAMtypetipo es estensione del file di output
Attenzione che la maggior parte dei programmi di allineamento sono creati con impostazioni di default basate sul genoma umano –> aggiustare i parametri sulla base del proprio organismo! Andiamo a lanciare il comando
# creiamo cartella dove salveremo file di allineamento
mkdir -p STAR/aligned
# lanciamo l'allineamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--genomeDir genome/index/ \
--readFilesIn Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz \
--readFilesCommand zcat \
--outFileNamePrefix STAR/aligned/SU_1_ \
--alignIntronMax 2000 \
--outSAMtype BAM SortedByCoordinateSe è tutto a posto l’allineamento dovrebbe impiegare circa 10 secondi. Come vedete la lunghezza massima media per gli introni di Arabidopsis è di circa 2000 bp, e gli ho detto che voglio che il file di output sia in formato .bam (formato allineamento binario). Il file .bam è in un formato che non può essere aperto con il terminale, ma che andremo poi a visualizzare con un apposito software. Per completezza vi inserisco un immagine dello stesso file di allineamento ma in un formato visualizzabile e consultabile sul terminale. Il formato di questo file è .sam (Sequence Alignment Format), formato non molto utilizzato perché non essendo compresso occupa molto spazio
SAM file
Andiamo ora a vedere i nostri file creati
# il file che ci interessa è chiamato SU_1_Aligned.sortedByCoord.out.bam
ls -rht STAR/aligned/
# leggiamo qualche statistica in questo altro file
less -S STAR/aligned/SU_1_Log.final.outQuanto vi risulta la percentuale di Uniquely mapped reads?
ALIGNMENT VISUALIZATION
Dopo aver effettuato l’allineamento può essere utile effettuare un controllo visivo sull reads allineate. A tal fine si utilizza spesso un software chiamato IGV (http://software.broadinstitute.org/software/igv/) che permette tra le altre cose di visualizzare i file dell’allineamento dopo aver caricato un genoma di riferimento
Integrative Genomics Viewer
Andiamo a scaricare il programma dal sito http://software.broadinstitute.org/software/igv/download cliccando in basso su “Launch with 2 GB”. Attendiamo lo scaricamento e apriamo il programma con doppio clic. Una volta aperto in alto a sinistra selezionare il genoma di A. thaliana TAIR10. A questo punto possiamo caricare i risultati del nostro allineamento, non prima però di aver creato un indice per il file .bam. Un altro indice direte voi, sì questa volta però non eseguiremo nessun passaggio ma vi farò prelevare il file dell’indice già preparato dal nostro server. Per completezza vi mostro il comando utilizzato per creare l’indice con Samtools
samtools index -b SU_1_Aligned.sortedByCoord.out.bamSu questi computer non possiamo installare il programma che si utilizza, tra le altre cose, anche per la creazione degli indici per i file di allineamento. Il programma si chiama Samtools (http://www.htslib.org/) ed è molto utilizzato in bioinformatica perché svolge numerose funzioni. Andiamo a copiare nella cartella dove sta il file .bam appena creato da noi, il file dell’indice, che permette che il file .bam venga visualizzato correttamente in IGV
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/
# il file che vedete deve avere lo stesso nome del file di allineamento, ma ha un'estensione diversa, .baiA questo punto potete tornare su IGV e cliccare su “File -> Load from file” e selezionare dalla cartella aligned il file SU_1_Aligned.sortedByCoord.out.bam. Selezioniamo Chr1 e zoomiamo sino a che la rotella inizia a girare, a quel punto possiamo visualizzare le nostre reads allineate ai geni del genoma di riferimento
Non sarà così facile trovare un po’ di reads!
Nel caso non funzionasse, copiare dal server alla vostra cartella, sia il file .bam che il .bai (potrebbe essere che il bai debba essere creato dallo stesso bam per funzionare la visualizzazione)
# bam
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam STAR/aligned/
# bai
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/La lunghezza degli introni è stata rispettata?
READS COUNTING
Una volta controllato che l’allineamento è andato a buon fine e siamo soddisfatti di come il programma di mapping ha lavorato, possiamo spostarci nel prossimo step verso l’analisi dei geni, ovvero la conta delle reads. Per conta delle reads si intende il contare quante reads sono state assegnate (hanno “mappato”) a un determinato gene (o trascritto, o esone, le cosiddette genomic features) in seguito al mapping. Per questo passaggio utilizzeremo il software featureCounts (http://subread.sourceforge.net/) del pacchetto Subread. Scarichiamo e testiamo il programma
# creiamo cartella
mkdir featureCounts
# scarichiamo dal server
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/featureCounts featureCounts/
# testiamo funzionamento
./featureCounts/featureCounts -vVediamo le opzioni e i parametri di featureCounts:
featureCountslanciamo il programma-aforniamo il file di annotazione del trascrittoma-opercorso e nome file di output-pindichiamo che lavoriamo su reads paired e che verranno contati i frammenti e non le reads-tindica quale feature contare-gindica che vogliamo raggrppare tutte le feature in questa meta-feature-Tquanti threads vogliamo usare per correre il programmain.bamnome del file di allineamento che vogliamo come input
A questo punto andiamo a lanciare il programma, ricordandoci che nel nostro caso vogliamo contare le feature del file di allineamento .bam
# creiamo cartella per output
mkdir -p featureCounts/output
# lanciamo la conta
./featureCounts/featureCounts \
-a genome/Arabidopsis_thaliana.TAIR10.39.gtf \
-o featureCounts/output/SU_1_counts.txt \
-p \
-t exon \
-g gene_id \
-T 10 \
STAR/aligned/SU_1_Aligned.sortedByCoord.out.bamDovrebbe terminare in qualche secondo. Andiamo ad esplorare i file di output
ls -lrht featureCounts/output/
less -S featureCounts/output/SU_1_counts.txt.summary
less -S featureCounts/output/SU_1_counts.txtLe colonne che ci interessano sono la prima e l’ultima (7)